home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byt1186b.arc / APPROX.LBR / SIMQ.C < prev   
Text File  |  1986-04-11  |  3KB  |  173 lines
  1. «PL85»
  2.  *
  3.  *    Solution of simultaneous linear equations AX = B
  4.  *    by Gaussian elimination with partial pivoting
  5.  *
  6.  *
  7.  *
  8.  * SYNOPSIS:
  9.  *
  10.  * double A[n*n], B[n], X[n];
  11.  * int n, flag;
  12.  * int IPS[];
  13.  * int simq();
  14.  *
  15.  * ercode = simq( A, B, X, n, flag, IPS );
  16.  *
  17.  *
  18.  *
  19.  * DESCRIPTION:
  20.  *
  21.  * B, X, IPS are vectors of length n.
  22.  * A is an n x n matrix (i.e., a vector of length n*n),
  23.  * stored row-wise: that is, A(i,j) = A[ij],
  24.  * where ij = i*n + j, which is the transpose of the normal
  25.  * column-wise storage.
  26.  *
  27.  * The contents of matrix A are destroyed.
  28.  *
  29.  * Set flag=0 to solve.
  30.  * Set flag=-1 to do a new back substitution for different B vector
  31.  * using the same A matrix previously reduced when flag=0.
  32.  *
  33.  * The routine returns nonzero on error; messages are printed.
  34.  *
  35.  *
  36.  * ACCURACY:
  37.  *
  38.  * Depends on the conditioning (range of eigenvalues) of matrix A.
  39.  *
  40.  *
  41.  * REFERENCE:
  42.  *
  43.  * Computer Solution of Linear Algebraic Systems,
  44.  * by George E. Forsythe and Cleve B. Moler; Prentice-Hall, 1967.
  45.  
  46.  
  47. int simq( A, B, X, n, flag, IPS )
  48. double A[], B[], X[];
  49. int n, flag;
  50. int IPS[];
  51. {
  52. int i, j, ij, ip, ipj, ipk, ipn;
  53. int idxpiv, iback;
  54. int k, kp, kp1, kpj, kpk, kpn;
  55. int nip, nkp, nm1;
  56. double em, q, rownrm, big, size, pivot, sum;
  57. double abs();
  58.  
  59. if( flag < 0 )
  60.     goto solve;
  61.  
  62. /*    Initialize IPS and X    */
  63.  
  64. ij=0;
  65. for( i=0; i<n; i++ )
  66.     {
  67.     IPS[i] = i;
  68.     rownrm = 0.0;
  69.     for( j=0; j<n; j++ )
  70.         {
  71.         q = abs( A[ij] );
  72.         if( rownrm < q )
  73.             rownrm = q;
  74.         ++ij;
  75.         }
  76.     if( rownrm == 0.0 )
  77.         {
  78.         puts("SIMQ ROWNRM=0");
  79.         return(1);
  80.         }
  81.     X[i] = 1.0/rownrm;
  82.     }
  83. /*                    
  84. /*    Gaussian elimination with partial pivoting     */
  85.  
  86. nm1 = n-1;
  87. for( k=0; k<nm1; k++ )
  88.     {
  89.     big= 0.0;
  90.     for( i=k; i<n; i++ )
  91.         {
  92.         ip = IPS[i];
  93.         ipk = n*ip + k;
  94.         size = abs( A[ipk] ) * X[ip];
  95.         if( size > big )
  96.             {
  97.             big = size;
  98.             idxpiv = i;
  99.             }
  100.         }
  101.  
  102.     if( big == 0.0 )
  103.         {
  104.         puts( "SIMQ BIG=0" );
  105.         return(2);
  106.         }
  107.     if( idxpiv != k )
  108.         {
  109.         j = IPS[k];
  110.         IPS[k] = IPS[idxpiv];
  111.         IPS[idxpiv] = j;
  112.         }
  113.     kp = IPS[k];
  114.     kpk = n*kp + k;
  115.     pivot = A[kpk];
  116.     kp1 = k+1;
  117.     for( i=kp1; i<n; i++ )
  118.         {
  119.         ip = IPS[i];
  120.         ipk = n*ip + k;
  121.         em = -A[ipk]/pivot;
  122.         A[ipk] = -em;
  123.         nip = n*ip;
  124.         nkp = n*kp;
  125.         for( j=kp1; j<n; j++ )
  126.             {
  127.             ipj = nip + j;
  128.             A[ipj] = A[ipj] + em * A[nkp + j];
  129.             }
  130.         }
  131.     }
  132. kpn = n * IPS[n-1] + n - 1;    /* last element of IPS[n] th row */
  133. if( A[kpn] == 0.0 )
  134.     {
  135.     puts( "SIMQ A[kpn]=0");
  136.     return(3);
  137.     }
  138. /*                
  139. /*    back substitution    */
  140.  
  141. solve:
  142. ip = IPS[0];
  143. X[0] = B[ip];
  144. for( i=1; i<n; i++ )
  145.     {
  146.     ip = IPS[i];
  147.     ipj = n * ip;
  148.     sum = 0.0;
  149.     for( j=0; j<i; j++ )
  150.         {
  151.         sum += A[ipj] * X[j];
  152.         ++ipj;
  153.         }
  154.     X[i] = B[ip] - sum;
  155.     }
  156.  
  157. ipn = n * IPS[n-1] + n - 1;
  158. X[n-1] = X[n-1]/A[ipn];
  159.  
  160. for( iback=1; iback<n; iback++ )
  161.     {
  162. /* i goes (n-1),...,1    */
  163.     i = nm1 - iback;
  164.     ip = IPS[i];
  165.     nip = n*ip;
  166.     sum = 0.0;
  167.     for( j=i+1; j<n; j++ )
  168.         sum += A[nip+j] * X[j];
  169.     X[i] = (X[i] - sum)/A[nip+i];
  170.     }
  171. return(0);
  172. }
  173.